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Recent developments in fly-by-wire control architectures for rotorcraft have introduced 
new interest in the identification of time- varying pilot control behavior in multi-axis control 
tasks. In this paper a maximum likelihood estimation method is used to estimate the 
parameters of a pilot model with time-dependent sigmoid functions to characterize time- 
varying human control behavior. An experiment was performed by 9 general aviation pilots 
who had to perform a simultaneous roll and pitch control task with time-varying aircraft 
dynamics. In 8 different conditions, the axis containing the time-varying dynamics and the 
growth factor of the dynamics were varied, allowing for an analysis of the performance of 
the estimation method when estimating time-dependent parameter functions. In addition, 
a detailed analysis of pilots’ adaptation to the time-varying aircraft dynamics in both the 
roll and pitch axes could be performed. Pilot control behavior in both axes was significantly 
affected by the time-varying aircraft dynamics in roll and pitch, and by the growth factor. 
The main effect was found in the axis that contained the time-varying dynamics. However, 
pilot control behavior also changed over time in the axis not containing the time-varying 
aircraft dynamics. This indicates that some cross coupling exists in the perception and 
control processes between the roll and pitch axes. 


Nomenclature 


Ad 

sinusoid amplitude, deg 

L 

likelihood function, - 

B 

bias 

M 

parameter function time of maximum growth, s 

Cqq 

Cramer-Rao lower bound 

Mqq 

Fisher information matrix 

e 

tracking error signal, deg 

m 

number of samples 

f 

probability density function, - 

N 

number of points 

fd 

disturbance forcing function, deg 

n 

pilot remnant signal, deg 

G 

parameter function growth rate, s^ 1 

n d 

forcing function frequency integer factor 

H c 

aircraft dynamics response 

Pi 

initial parameter function value, - 

H ol 

open-loop response 

Pi 

final parameter function value, - 

Hp 

pilot response 

s 

Laplace variable 

H n 

remnant filter 

T c 

vehicle dynamics time constant, s 

3 

imaginary unit 

Ti 

pilot lead time constant, s 

K c 

vehicle dynamics gain, - 

T m 

measurement time, s 

K n 

pilot remnant gain, - 

t 

time, s 

K v 

pilot visual gain, - 

u 

pilot control signal, deg 

k 

sinusoid index 
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Symbols 

a 

line search parameter, - 

a 2 * 

T e 

T v 

e 

prediction error, deg 

LO 

<t> 

roll angle, deg 

CO c 

<t>d 

sinusoid phase shift, rad 

nm 

tynn 

phase margin, deg 

Wd 

e 

parameter vector 

C nm 

e 

pitch angle, deg 


variance 

effective time delay, s 

pilot visual delay, s 

frequency, rad s _1 

crossover frequency, rad s^ 1 

pilot neuromuscular frequency, rad s _1 

sinusoid frequency, rad s -1 

pilot neuromuscular damping, - 


I. Introduction 


In many aircraft control tasks, pilots need to adapt their control strategy to controlled aircraft dynamics 
that change over time. Examples of control tasks with time-varying dynamics are: the manual control of a 
tilt-rotor aircraft during the transition from hovering to forward flight, and the manual control of rotorcraft 
with variable speed rotors or aircraft with reconfigurable flight control systems. 1 Currently, not much is 
known about the adaptation of pilot skill based control behavior to time- varying dynamics, as most previous 
research on control behavior focused on control tasks with constant dynamics. In addition, most previous 
research focused on pilot control behavior in single-axis control tasks, while most pilot control tasks require 
simultaneous control inputs in multiple axes. 

Recent developments in fly-by-wire control architectures for rotorcraft have introduced new interest in 
the identification of time-varying pilot control behavior in multi-axis control tasks. 2 Most research in the 
past on the time- varying aspect of human control behavior has been devoted to the development of complex 
mathematical models to describe adaptations of behavior over time. An accurate estimate of the parameters 
of these models would give valuable insight by quantifying how a human adapts to time-varying controlled 
dynamics. However, because of the complexity of the models it is not possible to accurately estimate their 
parameters from measured data. In addition, there are insufficient techniques for the identification of time- 
varying human control behavior and the estimation of time-varying pilot model parameters. 

Techniques for the identification and prediction of human control behavior in control tasks with constant 
dynamics are well established. 4-7 In recent years, work on the identification of time-varying pilot control 
behavior focused on the use of wavelet transforms to identify time- varying pilot frequency response functions 
and derived parameters such as crossover frequencies and phase margins. 8-10 However, when applied to 
human control data, the wavelet transform is known to produce inaccurate results, as it is very sensitive to 
the relatively high pilot remnant levels. 11 Furthermore, the use of wavelets does not allow for the estimation 
of time-varying pilot model parameters directly. A second step is required, fitting the pilot model frequency 
response to the wavelet frequency response. 6 

In a recent study, maximum likelihood estimation (MLE) has been used to estimate time-varying pilot 
model parameters directly from time-domain data. 11 The time dependency of the pilot model parameters was 
introduced by estimating constant parameters on a sliding time window. This approach produced relatively 
noisy results and required significant computational power as the calculations have to be performed at every 
time step. 

In this paper an MLE method is used to estimate the parameters of a pilot model with time-dependent 
sigmoid functions to characterize time-varying human control behavior from experiment data in a single 
step. This method promises to be far more efficient to the time-windowing MLE method previously used 
and possibly allows for more accurate results. An experiment was conducted using 9 general aviation pilots 
who performed a simultaneous roll and pitch control task with time- varying aircraft dynamics. In 8 different 
conditions the axis containing the time-varying dynamics and growth factor of the dynamics were varied, 
allowing for an analysis of the performance of the MLE method when estimating time-dependent parameter 
functions. In addition, a detailed analysis of pilots’ adaptation in different axes of control to the time- varying 
aircraft dynamics was performed. 

The multi-axis control task and the time-varying aircraft and pilot model dynamics will be discussed in 
Sec. II, followed by the adopted MLE method used for the estimation of time- varying pilot model parameters 
in Sec. III. Next, the experiment setup will be discussed in Sec. IV. The results will be given in Sec. V and 
discussed in Sec. VI. Finally, conclusions will be given in Sec. VII. 
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II. Multi-Axis Control Task 


A. Control Loop 

Fig. 1 depicts the multi-axis disturbance-rejection task that was used to analyze time-varying pilot control 
behavior in this study. A pilot controls the roll and pitch attitude (</> and 6) of an aircraft with dynamics H c 
by making control inputs ( u ^ and ug) with a joystick. The aircraft dynamics do not have any cross coupling 
between the roll and pitch axes. The pilot’s task is to keep the aircraft level by keeping the roll and pitch 
angles (c t> and 9) as small as possible at all times, while they are perturbed by disturbances (/d^ and fgg)- 
The disturbed roll and pitch angels (</>d and Og) are displayed on the primary flight display depicted in Fig. 2. 

In single-loop control tasks, pilot skill-based control behavior can be modeled by a linear response function 
H p and a remnant signal n that accounts for non-linear behavior. As the aircraft dynamics in roll and pitch 
do not have any cross coupling, it is assumed that pilot control behavior in this task can be modeled by two 
independent linear response functions {H pt f, and H p g) and two independent remnant signals (rig, and rig). 
The resulting control diagram structure of the pilot-vehicle system consists of two independent single-loop 
control loops (Fig. 1). Hence, the system can be termed multiple single-loop. The assumption of two separate 
control-loops in a control task with uncoupled dynamics has been adopted in previous studies and was found 
to be appropriate for data analysis purposes. 12 



Figure 1. Multi-axis closed-loop control task. 



A typical example of the multi-axis control task depicted in Fig. 1 is longitudinal and lateral stabilization 
of an aircraft in straight, wings-level, horizontal flight using aileron and elevator control inputs. Techniques 
for the identification and parameterization of pilot response functions in this type of control task are well 
established for vehicle dynamics H c that are constant over time. 4, 13,14 However, the vehicle dynamics in 
Fig. 1 change over time. Such time-varying dynamics can, for example, be introduced by an aircraft control 
system that adapts to different flight conditions or phases of flight. As a result, pilots will adapt their control 
behavior over time to maintain adequate levels of performance and stability margins. The goal of this study 
is to quantify this adaptation in pilot control behavior by estimating time-varying pilot model parameters 
using maximum likelihood estimation. 11 
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B. Time-Dependent Parameter Functions 

Many types of parameter functions can be used to describe a time dependency in the pilot model parameters. 
However, the function needs to be able to describe the expected transient response of a certain parameter. 
The most appropriate parameter function needs to be determined with prior knowledge of the change in 
controlled dynamics and by carefully inspecting the measured data. In practice, the complexity of the func- 
tion will depend on many variables, such as the overshoot and settling time of the transient response. Care 
should be taken that not too many additional parameters are introduced by the parameter function, as this 
can lead to over parameterization of the complete pilot model, resulting in less accurate parameter esti- 
mates. Furthermore, it will be more difficult to find a global optimum parameter set, as multiple parameter 
combinations may result in similar model responses. 

The parameter function considered in this study is the generalized logistic function, a widely-used sigmoid 
function. Sigmoid functions are characterized by their S-shape, describing the transition from an initial value 
to a final value. These functions are often used for learning curve and growth modeling. The generalized 
logistic function as a function of time t is defined by: 


P(t) =Pi + 


P 2 -P 1 

1 + g-G(t-M) 


(1) 


where Pi is the initial pilot model parameter value and P 2 the final parameter value. The time of maximum 
growth is defined by M and the growth rate by G. Fig. 3 depicts the generalized logistic function for 
variations of the time-of-maximum-growth and growth-rate parameters. For G —> 00 the parameter function 
converges to a step function, allowing for the modeling of an instant change in pilot control behavior. Note 
that for a value of G = 10 the function already approaches a step function (Fig. 3b). This means that this 
parameter might be difficult to estimate, depending on its real value and the sampling frequency of the data, 
as it has no upper limit (e.g. a growth rate of 100 or 1000 are indistinguishable). 


(a) M variation (Pi = 0.1, P 2 = 0.4, G = 0.5) 


(b) G variation (Pi — 0.1, P 2 = 0.4, M = 50) 




Figure 3. Generalized logistic function parameter variation. 


C. Time-varying Dynamics 

The time- varying aircraft roll and pitch dynamics in Fig. 1 are defined by: 


H c (s,t) 


K c {t) 

s 2 + T c (t)s 


( 2 ) 


where the gain K c (t) = K c \ + [if c2 — K c i]/[l + e - G ( ( - M )] and the time constant T c (t) = T c i + [T c 2 — T c i]/[1 + 
e -G(t-M)] are time-varying parameters defined by the parameter function given in Eq. (1). The initial and 
final parameter function values are K c i = 90.0, K C 2 = 30.0, T c i = 6.0 s, and T c2 = 0.2 s. The time of 
maximum growth is M = 50.0 s and the growth rate is either G = 0.5 s^ 1 or G = 100.0 s -1 , as will be 
explained in Sec. IV. The time of maximum growth and the growth rate are equivalent for both time- varying 
parameters of the aircraft dynamics. Given the parameter functions described above, the aircraft dynamics 
will change over time from H c i = 90.0/(s 2 + 6.0s) to H c 2 = 30.0/(s 2 + 0.2s). 
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(a) Aircraft dynamics magnitude 


(b) Aircraft dynamics phase 




lcr 1 io° 10 1 


ui, rad s 1 


Figure 4. Frequency response of the aircraft dynamics at different time instances. 


Fig. 4 depicts the time-varying aircraft dynamics for different time instances. Note that the dynamics 
will change from mostly stable dynamics (1/s) to mostly unstable dynamics (1/s 2 ) in the human crossover 
frequency range (1-5 rad/s). The crossover model theorem states that human operators adapt their control 
behavior such that the pilot-vehicle open-loop response in the crossover region can be described by: 

H 0 l(s) = H p (s)H c (s) « ye- ST * (3) 

with uj c the crossover frequency and r e a time delay lumping the human processing lags. 4 In order to comply 
with Eq. (3) and given the variations in H c , the equalization dynamics of the pilot model need to contain 
a lead time constant. The pilot lead time constant is expected to increase over time (lead is generated at 
increasingly lower frequencies) as the time constant of the aircraft dynamics decreases. Based on previous 
research, the pilot gain is also expected to change over time. 10 The time- varying pilot model capable of 
modeling this adaptation in pilot control behavior is defined by: 


H p {s,t ) 


equalization 


limitations 

✓n. 


K v {t)[l + Ti(t)s] e- 




2C 


nm^nm 


( 4 ) 


The equalization dynamics of the pilot model consist of a time- varying gain K v (t ) = K v \ + [K v 2 — 
K v i]/[l + e~ G ( t-M )] and a time- varying lead time constant Tj(f) = Tn + [7] 2 — Tii]/[l + e~ G ^~ M ^]. The pilot 
limitations, defined by the time delay t v , neuromuscular damping £ nrn , and neuromuscular frequency u> nm 
- are assumed to be constant. Previous research found that the limitation parameters are approximately 
equivalent for a pilot controlling dynamics similar to H c \ or similar to H c 2 . 15 The time of maximum growth 
and the growth rate of both parameter functions are assumed to be equal. This means a parameter vector 
0 = [K v i K v2 Tn T 12 M G t v ( nm ui nm ] with a total of 9 parameters needs to be estimated by the parameter 
estimation procedure. 


D. Disturbance Forcing Functions 

Both disturbance forcing functions are a sum of sines constructed using the following equation: 


N d 

fdcj>,do(t) = ^2 A d(k) sin[uj d (k)t + <t>d 4 >,de(k)] (5) 

k = 1 

with Nd = 10 the number of sine waves, and u)d, Ad and 4>d<i>,de the frequency, amplitude and phase shift 
of the A: th sine wave, respectively. The length of an experiment run is T = 90.00 s. The measurement time 
used to construct the forcing functions is T m = 81.92 s. With a data sampling frequency of 100 Hz, this 
measurement time contains the highest power-of-two measurements (N m = 8192) in the total length of an 
experiment run. The sinusoid frequencies tOd(k) were all integer multiples of the measurement-time base 
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Table 1. Disturbance forcing function properties. 


fc, - 

rid, ~ 

LJd, rad s 1 

A d , deg 

rad 

4>d.e, rad 

1 

3 

0.230 

1.186 

0.974 

-0.753 

2 

5 

0.384 

1.121 

-3.026 

1.564 

3 

8 

0.614 

0.991 

0.744 

0.588 

4 

13 

0.997 

0.756 

-2.300 

-0.546 

5 

22 

1.687 

0.447 

2.984 

0.674 

6 

34 

2.608 

0.245 

-2.513 

-1.724 

7 

53 

4.065 

0.123 

2.211 

-1.963 

8 

86 

6.596 

0.061 

1.004 

-2.189 

9 

139 

10.661 

0.036 

-2.255 

0.875 

10 

229 

17.564 

0.025 

2.210 

0.604 


frequency, w m = 27r/T m = 0.0767 rad/s, and were covering the frequency range of human control (0.1-18 
rad/s). The spacing between the frequencies was approximately equal on a logarithmic scale. 

The absolute value of a second-order low-pass filter at a sinusoid frequency was used to determine the 
amplitudes of the individual sine waves: 


A i{k) 


[1 + Q.ljujdjk)} 2 
[1 + 0.8 juj d {k)f 


( 6 ) 


The second-order low-pass filter ensures reduced amplitudes at higher frequencies, yielding a control 
task that is not overly difficult. The final amplitude distributions were scaled to produce roll and pitch 
disturbance forcing functions with a standard deviation of 1.5 deg. 

To determine the forcing function phase distributions, a large number of random phase sets were gener- 
ated. The set yielding a signal with a probability distribution closest to a Gaussian distribution, without 
leading to excessive peaks, was selected. 16 The characteristics of the forcing functions are summarized in 
Table 1. Note that the the roll and pitch disturbance forcing functions contain the same frequencies and 
amplitudes. However, the phase distributions of the two forcing functions were different. 


III. Identification of Time- Varying Pilot Control Behavior 

A. Maximum Likelihood Estimation 

The MLE procedure adopted in this research is similar to the one described in Ref. 7. Some modifications 
are applied to facilitate the time-dependent parameter functions of the pilot model. 

The MLE procedure calculates an estimate 0 of the pilot model parameter vector 0 by maximizing 
a likelihood function. The likelihood function L (0) is defined as the joint conditional probability density 
function of the prediction error for in measurements of the pilot control signal u{k)\ 


L (0) = /[e(l),e(2),...,e(fc),...,e(m)|0] (7) 

The prediction error, indicated as e(fc) in Eq. (7), is defined as the difference between the measured pilot 
control signal u(k) and the control signal from the pilot model u(k ) at discrete instants. When the remnant 
is assumed to be an additive zero-mean Gaussian white noise signal, the conditional probability density 
function for one measurement of e(fc) is given by: 


f«k) |0) 


1 


e 2 (fc) 

2<x2 


( 8 ) 


The set of parameters that maximizes the likelihood function is the maximum likelihood estimate of 
the parameter vector 0. To create a more straightforward optimization problem, it is common practice 
to minimize the negative natural logarithm of the likelihood function instead of maximizing L(0). The 
parameter vector resulting in the global minimum of the negative log-likelihood function is then the maximum 
likelihood estimate ©ml • For the single-output pilot model defined in Eq. (4), combining Eq. (7) and Eq. (8) 
results into the following expression for the maximum likelihood parameter estimate: 
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Oml — argmin —In L (0) = argmin 
e e 


m. 2 

¥ ln 


m 

M^ c2(fc) 


k = 1 


( 9 ) 


The minimum of this function is found using an unconstrained gradient-based Gauss-Newton optimiza- 
tion. The iterative parameter update equation for the Gauss-Newton optimization is given by: 


0 (* + 1) = 0 (i) - a (i) Mee [© (i)] (10) 

where a is the line-search parameter determined at each iteration before the actual parameter update. The 
optimal a is determined by calculating the likelihood for different values of a between 0 and 1. The a that 
results in the biggest decrease of the negative natural logarithm of the likelihood function is used for the 
parameter update. This ensures the most rapid minimization of the likelihood function. 

The gradients of the likelihood function with respect to all parameters, dL/dO, are calculated using a 
simple three-point finite difference approximation: 


dL[Q (i)] _ LjO (i) + A0] - L[Q (i) - A0] 

do 2A0 1 ’ 

where A0 = 0.01 is a small increment to the parameters in the parameter vector. The Fisher information 
matrix, indicated with the symbol Mqq in Eq. (10), is given by: 


7 11 

Mee = Y, 

n fc = i 

The Fisher information matrix is symmetrical and needs to be positive definite (i.e., full rank), as it 
needs to be inverted for the parameter update equation. The inverse of the Fisher information matrix yields 
the Cramer-Rao lower bound (CRLB), which determines the lowest achievable variance of the parameter 
estimates: 


de(k) 

do 


(12) 


Gee = Mqq (13) 

Both the Fisher information matrix and the first-order gradient of the likelihood function vary with the 
magnitude of the estimation error e, yielding parameter update steps of larger magnitude if estimation errors 
are larger. 

As described in Ref. 7, a genetic algorithm is used to find a solution of the parameter vector close to the 
global minimum. This solution is then used as an initial parameter vector in the unconstrained gradient- 
based Gauss-Newton optimization. This significantly increases the chance of finding the global optimum 
solution of the parameter estimation problem. 

B. Bias and Variance 

A Monte Carlo analysis with 140 simulations was performed to determine the bias and variance characteristics 
of the proposed MLE method in combination with the control task and aircraft dynamics discussed in Sec. II. 
The simulations were performed using a single control loop of the diagram in Fig. 1. Note that both the roll 
and pitch control loops in the diagram are equivalent. Fig. 5 depicts the simplified control loop used in the 
simulations. 

The time-varying aircraft dynamics H c and pilot dynamics H p used in the simulations are described in 
Sec. II. C. The growth rate of the time-varying aircraft dynamics was set to 0.5. The parameters of the 
time-varying pilot model were selected using data from a test experiment: K v i = 0.09, K v 2 = 0.07, Tn = 0.4 
s, T 12 = 1, 2 s , M = 50.0 s, G = 0.5 s -1 , r v = 0.28 s, ( nm = 0.35, and w nm = 11.25 rad/s. 

The pilot remnant ng was simulated using a low-pass filtered white noise signal w. The remnant filter 
was defined by: H n = K n /(0.2s + 1). Three different levels of pilot remnant were used to investigate the 
effect of remnant intensity on the bias and variance of the parameter estimates. The remnant intensity is 
defined by: P n = cr 2 (n)/cr 2 (u). The filter gain K n was selected to induce remnant intensities of P„ = 0.05, 
P n = 0.15, and P n = 0.25. 

Figs. 6a and 6b depict the relative bias and variance, respectively, of the pilot model parameters for 
different levels of remnant intensity. The relative bias and variance are calculated by dividing the bias and 
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Figure 5. Single-loop control task used for simulations. 


variance of each parameter by the value of the parameter itself. From Fig. 6a it can be observed that the 
bias increases for most of the parameters as the remnant intensity increases. Furthermore, the bias of K v 2 , 
T 12 , and G is higher compared to the other variables. 

Fig. 6b indicates that the relative variance of the parameters also increases as the remnant intensity 
increases, an expected result. The figure further indicates that the variance of the final parameter function 
values K V 2 and T 12 is higher compared to the initial values K v 1 and Tn. Finally, it can be observed that the 
growth rate G has the highest variance of all parameters. 

Table 2 provides the CRLB variance of the pilot model parameters averaged over 140 simulations for 
different values of remnant intensity. Although the values in this table provide the lowest achievable variance, 
they still provide a good indication of the relative estimation accuracy of every pilot model parameter. 

Table 2. CRLB variance of the time-varying pilot model parameters (140 simulations). 


Ce© 


Pn 

K v i 

K v 2 

T n 

Tn 

M 

G 

T v 

Cnm 

UJnm 

0.05 

5.49 Xl0 -6 

1.56xim 5 

3.78X10 -4 

4.67 X 10 -3 

3.48X10 -1 

3.62X10 -2 

1.71X10 -5 

2.79X10 -4 

7.02X10 -2 

0.15 

4.81xl0~ 6 

1.27X10 -5 

3.49X10 -4 

4.99X10 -3 

3.67X10 -1 

3.90X10+ 6 

1.56x 10 -5 

2.57 X 10 -4 

6.87X10 -2 

0.25 

3.77xl0~ 6 

l.lOx 10~ 5 

2.85X10” 4 

8.52X10" 3 

5.96X10 -1 

4.86X10+ 7 

1.31X10 -5 

2.07 X 10~ 4 

6.08X 10 -2 


The table indicates that the final parameter function values have higher lower bounds compared to the 
initial parameter function values, implicating that these variables possibly have a lower accuracy. This was 
also apparent from Fig. 6b. Next, the growth rate appears to have the highest lower bound values of all 
model parameters. The low expected accuracy of this parameter was also apparent from Fig. 6b. 

Summarizing, the results of the Monte Carlo analysis indicate that the final parameter function variables 
K v 2 and Tj 2 , and the growth rate G will most likely have the lowest accuracy when estimated from experiment 
data. However, the results presented are for this particular set of pilot model parameters. The true pilot 
model parameter bias and variance characteristics will depend on many factors, such as pilot performance, 
remnant characteristics, and control behavior. 


(a) Pilot model parameter bias 


(b) Pilot model parameter variance 



parameter 



K v 1 K v2 T n T12 M G 
parameter 


Figure 6. Relative bias and variance of the time- varying pilot model parameters (140 simulations). 
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IV. Experiment Setup 


The experiment was performed with two goals in mind; to investigate how pilot control behavior is 
affected by time-varying aircraft dynamics in a multi-axis control task and to evaluate the performance of 
the proposed maximum likelihood estimation method when estimating the parameters of a time- varying pilot 
model with time-dependent parameter functions in different controlled axes and for different growth rates of 
the controlled aircraft dynamics. 

A. Method 

1. Independent Variables 

Three independent variables with two levels each were varied in the experiment: 1) roll-axis dynamics (H c i 
or H c i — > H c 2), 2) pitch-axis dynamics (H c 1 or H c 1 — > H c 2), and 3) growth rate of the change in dynamics 
( G = 0.5 or G = 100.0). A full- factorial design was implemented, resulting in 8 conditions in total. However, 
two of these conditions are equivalent, as the growth rate does not have an effect when the dynamics in both 
the roll and pitch axes remain constant {H c 1 ). This means only 7 distinctive conditions result from this 
full-factorial design. As a reference, an additional condition was performed with constant dynamics H c 2 in 
both the roll and pitch axes. The aircraft dynamics H c 1 and H c 2 are defined in Sec. II. C. 

The 8 experiment conditions are listed in Table 3. Conditions 1 and 2 - referred to by the symbols Cl and 
C2 - have constant dynamics in both the roll and pitch axes. Conditions 3 to 5 (C3-C5) have time-varying 
dynamics with a growth rate of 0.5 in one or both of the axes. This growth rate results in a gradual change 
from H c i to H c 2 in about 20 s. In conditions 6 to 8 (C6-C8), the axis containing the time- varying dynamics 
are equivalent to C3 to C5. However, the time-varying dynamics have a growth rate of 100.0, resulting in 
an almost instant (step- like) change from H c 1 to H c i. 


Table 3. Experiment conditions. 


Condition 

Roll dynamics 

Pitch dynamics 

Growth rate 

Cl 

H c i 


Hci 



C2 

h c2 


H c2 



C3 

H c i - 

+ -Hc2 

Hci 


0.5 X 

C4 

H C 1 


Hci - 

H c 2 

0.5 X 

C5 

H c i - 

+ Hc2 

Hci - 

H c 2 

0.5 X 

C6 

H c i - 

+ -Hc2 

Hci 


100.0 1 

C7 

H c i 


Hci - 

H c 2 

100.0 1 

C8 

-Hci - 

+ -Hc2 

He i - 

+ H c 2 

100.0 1 


2. Apparatus 

Fig. 2 shows the primary flight display (PFD) used to provide pilots with attitude information. The PFD 
was depicted on a ViewSonic G225f CRT monitor. This monitor has a viewing area with a diagonal of 20 
inch and was set to a resolution of 1600x1200 px. The size of the PFD was 400x400 px and was depicted in 
the center of the screen. The thickness of the aircraft symbol was 6 px and the thickness of the horizon was 
2 px. The roll angle was magnified by a factor of 2 on the display. One degree of pitch angle was visible by 
an 8 px translation of the horizon from the aircraft symbol (see Fig. 2). Pilots were seated approximately 
1.8 ft in front of the monitor. 

A Logitech Extreme 3D Pro joystick was used to make control inputs. This right-handed joystick has con- 
siderable breakout force. Joystick output values in the roll and pitch axes ranged from -1 to 1, corresponding 
to full left-right and fore-aft joystick deflections, respectively. 

3. Participants and Instructions 

Nine general aviation pilots participated in the experiment. All were male with an average age of 24.3 years. 
The average total of flying hours was 773.2 hours with an average of 78.2 flying hours in the last 6 months. 
All pilots were right handed and only 1 had previous experience with the type of control task used. 

Every pilot received a briefing before the start of the experiment. The briefing contained information 
about the objective of the study. Furthermore, they were given details about the control task, conditions, 
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and procedures of the experiment. Pilots were instructed to make control inputs in an active and continuous 
manner, trying to minimize the roll and pitch angles presented on the visual display to the best of their 
capabilities. In addition, they were told to focus on minimizing roll and pitch angles at the same time. 

4- Procedures 

All experimental runs had a length of 90 s. The first 8.08 s of the run were considered the run-in time. 
During this phase, pilots were able to stabilize the controlled dynamics and adjust to the task. Data from 
the last 81.92 s were used as measurement data for the data analysis (see Sec. II. D). 

All pilots performed a total of 96 runs. The experiment consisted of 4 segments, each containing 24 
runs. Each of the 8 conditions was presented 3 times in every segment. The conditions were presented in a 
quasi-random order using a Latin-square design. Pilots were not informed about which condition they were 
performing. 

The first two segments were used for pilot training and runs from the last two segments were used for 
data analysis. After every run, tracking performance - defined as the root mean square (RMS) of the roll 
and pitch angles - was presented to the pilots to motivate them to constantly perform at their maximum 
level of performance. In between each segment, the experimenter made sure pilots reached and maintained 
asymptotic performance; that is, more or less constant RMS values. 

Pilots were able to go through the experiment at there own pace. They were able to start the next run 
by pressing the trigger on the joystick. This allowed them to take breaks when feeling any discomfort. Every 
pilot was able to complete the experiment within 4 hours. 

5. Dependent Measures 

Several dependent measures were calculated from the measured signals from the experiment. First, pilot 
performance and control activity in the roll and pitch axes of the disturbance-rejection task were evaluated 
in terms of the RMS of the attitude and control signals, respectively. Next, the parameters of the time- 
varying pilot model were estimated using the MLE method discussed in Sec. III. For every participant and 
condition the experimental data were averaged over six runs before the start of the estimation procedure. 
The control signal variance accounted for (VAF) was calculated as a measure for the accuracy of the pilot 
model in describing the measured control signal data. Finally, time- varying pilot-vehicle open-loop responses 
and corresponding time-varying crossover frequencies and phase margins were calculated. 

B. Hypotheses 

Very little research has been performed on the identification of time-varying control behavior in multi- 
axis control tasks, or the estimation of time-varying parameters using experiment data. However, some 
hypotheses on the performance of MLE for the estimation of time- varying parameters from experiment data 
can be formulated based on Sec. III.B. In addition, some hypotheses on pilot adaptation to the time-varying 
aircraft dynamics can be formulated based on experiments with single-loop control tasks with constant 
dynamics. 15, 17 

In Sec. III.B it was found that the bias and variance of the final parameter function parameters and the 
growth rate parameter were higher compared to the other parameters of the model. The bias and variance 
increased along with levels of remnant. It is expected that in a multi-axis task, pilot remnant levels will 
be higher, as attention needs to be divided between the two axes. This results in a hypothesis that larger 
variations are expected to be found in K v % , Tl2, and G. Furthermore, it is hypothesized that the accuracy 
of the overall pilot model (determined by the VAF) will be higher for the unstable dynamics H c 2 . 17 

In conditions in which the aircraft dynamics change from stable dynamics H c 1 to unstable dynamics H c 2 , 
it is hypothesized that pilot roll and pitch disturbance- rejection performance will decrease, and roll and pitch 
control activity will increase over time. 15 In addition, crossover frequencies will increase and phase margins 
will decrease. Finally, the pilot visual gain K v and visual lead time constant Tj are expected to decrease and 
increase, respectively, when the dynamics change from H c i to H c 2 . 15 
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(a) Roll attitude signal 


(b) Roll control signal 




condition condition 

Figure 7. Roll tracking performance and control activity. 


V. Results 

This section provides the combined results of the 9 pilots that participated in the experiment. The error 
bars in the error bar plots are corrected for between-subject variability by normalizing the subject means 
across conditions. 18 The data from the reference conditions Cl and C2 are marked by gray continuous and 
dashed lines, respectively. All the dependent measures of the experiment were analyzed using a repeated- 
measures analysis of variance (ANOVA) to uncover any significant trends. 

A. Pilot Performance and Control Activity 

Pilot performance and control activity in the roll axis - defined by the RMS of the roll angle and the RMS 
of the roll control signal - are depicted in Figs. 7a and 7b. Figs. 8a and 8b provide the same data for the 
pitch axis. A lower value for the roll or pitch angle RMS indicates a higher performance and a lower value 
for the control signal RMS indicates a lower control activity. 

Fig. 7a indicates that performance is higher for the reference condition Cl with stable dynamics compared 
to reference condition C2 with unstable dynamics, a highly significant effect [P(l, 8) = 49.661, p < 0.001]. 
This was an expected result, also observed in previous studies. 15 When looking at the effects of the change in 
roll dynamics on the RMS of the roll angle, a significant deterioration in performance can be observed when 
there is a change in roll dynamics from H c i to H c 2 [P(l, 8) = 41.302, p < 0.001]. The performance in these 
conditions was approximately half way between the performance for Cl and C2. No significant effect was 
found from a change in pitch dynamics [F(l, 8) = 3.252, p > 0.05]. The change in growth factor introduced 
a significant effect [F(l, 8) = 10.757, p < 0.05]. For the higher growth factor of the aircraft dynamics, the 
performance was slightly worse compared to the lower growth factor. One significant interaction was present 
between roll dynamics and growth factor [F(l, 8) = 15.441, p < 0.05]. For the higher growth factor, the 
change in roll dynamics introduced a larger decrease in roll performance. 

From Fig. 7b the observation can be made that the control activity is significantly higher for condition 
C2 compared to Cl [F(l, 8) = 39.329, p < 0.001]. A change in roll dynamics significantly increased roll 
control activity [_F(1, 8) = 73.000, p < 0.001]. Although roll control activity also seems slightly increased 
when only pitch dynamics vary in time (C4 and C7), the overall effect is not significant [F( 1, 8) = 2.189, 
p > 0.05]. Growth rate did not introduce a significant effect [F(l, 8) = 4.419, p > 0.05]. Two interactions 
between the independent variables were found. The first, between growth rate and roll dynamics [P(l, 8) 
= 16.115, p < 0.001]. This indicates that for the higher growth rate, the control activity increase from the 
change in roll dynamics is larger. The second interaction occurred between the change in roll and the change 
in pitch dynamics, reflecting the previously observed increase in control activity when only pitch dynamics 
vary in time [F(l, 8) = 9.564, p < 0.05]. 

The results in the pitch axis when the pitch dynamics were changed from stable to unstable dynamics are 
very similar to the results in the roll axis for a change in roll dynamics. From Fig. 8a it can be observed that 
disturbance-rejection performance is significantly better in the condition with constant stable dynamics Cl 
compared to the condition with constant unstable dynamics C2 [F(l, 8) = 29.778, p < 0.01]. The changing 
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(a) Pitch attitude signal 


(b) Pitch control signal 




condition condition 

Figure 8. Pitch tracking performance and control activity. 


dynamics in pitch [(C4,C5) and (C7,C8)] introduced a significant effect in pitch RMS [F(l, 8) = 29.004, 
p < 0.01], but no significant effect was found from the changing roll dynamics [F(l, 8) = 1.566, p > 0.05]. 
The growth rate of the changing aircraft dynamics also affected the pitch disturbance-rejection performance 
significantly [F(l, 8) = 6.406, p < 0.05]. The performance is slightly lower for the dynamics with a higher 
growth rate. One significant interaction is present in the data between growth rate and pitch dynamics [F(l, 
8) = 17.286, p < 0.01]. The reduction in performance caused by the change in pitch dynamics is slightly 
larger for the dynamics with a higher growth rate. 

Fig. 8b indicates that pitch control activity is significantly higher for the condition with constant unstable 
dynamics compared to the condition with constant stable dynamics [F(l, 8) = 29.882, p < 0.01]. The change 
in roll dynamics did not introduce a significant effect [F(l, 8) = 0.463, p > 0.05]. However, the change in pitch 
dynamics influenced the pitch control activity significantly [F(l, 8) = 36.488, p < 0.001]. In the conditions 
with time-varying dynamics in pitch, the control activity in pitch was significantly higher, approximately 
halfway between the control activity for Cl and C2. In addition, pitch control activity was significantly 
higher in the conditions with the higher growth rate [F(l, 8) = 17.670, p < 0.01]. The larger increase in 
control activity for the higher growth rate when the pitch dynamics changed over time was revealed by a 
significant interaction between growth rate and pitch dynamics [F(l, 8) = 11.843, p < 0.01]. 

Finally, by comparing performance and control activity between the roll and pitch axes (comparing Figs. 
7 and 8), it can be observed that the performance in reducing the disturbance was higher and control activity 
was lower in the pitch axis. These are both significant effects [F(l, 8) = 22.218, p < 0.01] and [F(l, 8) = 
50.480, p < 0.001]. 

B. Pilot Model Variance Accounted For 

The pilot model parameters were estimated using the MLE method discussed in Sec. III. For the con- 
ditions with constant dynamics (Cl and C2), a pilot model was identified without any time-dependent 
parameter functions. This simplified the model and reduced the parameter vector to 5 parameters, 0 = 
[K v i Tn T~v Cnm ^ nm ]■ 

Figs. 9 and 10 give the VAF for the roll and pitch control signal, respectively. A higher VAF indicates 
that a larger percentage of the measured data can be explained by the linear pilot model. In addition to 
the model with time-dependent parameter functions [H p (s,t)], the VAF in the conditions with time- varying 
dynamics was also calculated for the model with constant parameters [i7 p (s)]. This allows to determine if 
the pilot model with time-dependent parameter functions is able to describe the data better compared to 
the pilot model without parameter functions when the aircraft dynamics change over time. 

Fig. 9 reveals, as hypothesized in Sec. IV. B, that the VAF in the roll axis is significantly higher for the 
condition with unstable dynamics C2 compared to the condition with stable dynamics Cl [F(l, 8) = 25.590, 
p < 0.01]. An investigation of the effects of the change in roll and pitch dynamics and growth rate shows that 
only a change in roll dynamics introduced a significant effect in the roll VAF. The roll VAF in the conditions 
with time- varying roll dynamics was significantly higher [F(l, 8) = 80.525, p < 0.001], approximately at the 
same level of the VAF in condition C2. Comparing the VAF of both pilot models in the conditions with 


12 of 20 


American Institute of Aeronautics and Astronautics 




condition 


condition 


Figure 9. Roll axis variance accounted for. 


Figure 10. Pitch axis variance accounted for. 


time- varying aircraft dynamics (white compared to black squares in conditions C3-C5 and C6-C8), it can be 
observed that the VAF is slightly higher for the pilot model with time-dependent parameter functions [F(l, 
8) = 51.867, p < 0.001]. This indicates that this model is indeed more suitable to describe the data in the 
conditions with time- varying aircraft dynamics. 

The results in the pitch axis (Fig. 10) are very similar compared to the roll axis. The VAF in C2 is 
significantly higher compared to Cl [F(l, 8) = 18.329, p < 0.01]. The change in pitch dynamics from stable 
to unstable dynamics introduced a significant increase in VAF [F(l, 8) = 40.775, p < 0.001]. No significant 
effects were introduced by the change in roll dynamics or the growth rate. When comparing the two pilot 
models in the conditions with time- varying aircraft dynamics, it can be observed that the VAF is significantly 
higher for the model with time-dependent parameter functions [F(l, 8) = 91.753, p < 0.001]. Again, this 
indicates that this model is more suitable to describe the data in the conditions with time-varying aircraft 
dynamics. 

Finally, comparing the VAF between the roll and pitch axes (Fig. 9 compared to Fig. 10), the observation 
can be made that the VAF in the pitch axis was significantly higher [F(l, 8) = 107.690, p < 0.001]. This 
might indicate that pilots behaved more linear (i.e., introduced less remnant) in the pitch axis, as a higher 
percentage of the data could be explained by the linear pilot model. 

C. Pilot Model Parameters 

Figs. 11 and 12 present the estimated parameters for the pilot model in roll and pitch, respectively. The 
parameters estimated in all conditions are depicted by black squares. The parameters that are only available 
in the conditions with time- varying aircraft dynamics are depicted by white squares. As for these parameters 
a full-factorial repeated- measures AN OVA cannot be performed, a repeated-measures AN OVA is performed 
with the time-varying aircraft dynamics (not taking into account roll and pitch) and growth rate as factors. 

Results from the ANOVA indicate that, for the roll axis, there is no significant difference in the initial 
visual gain K v \ between the condition with constant stable dynamics Cl and the condition with constant 
unstable dynamics C2 [F(l, 8) = 2.159, p > 0.05] (see Fig. 11a). No significant effects are found for K v \ for 
the change in roll or pitch dynamics, or the growth rate of the dynamics. This is an expected result, as K v \ 
characterizes pilot control behavior for the stable aircraft dynamics H c \ in all conditions, except in C2. For 
the final visual gain I\ v 2 a significant effect was found for the change in dynamics [F( 2, 16) = 5.889, p < 
0.05]. This was also an expected result, because the values for this parameter in conditions with a change in 
roll dynamics (C3,C5,C6,C8) should be closer to the value for K vi in C2. Comparing K vl with K v 2 for the 
conditions with time-varying dynamics, a significant lower value for K v 2 was found [F(l, 8) = 8.604, p < 
0.05]. 

The visual gain for the pitch pilot model is given in Fig. 12a. The initial visual gain K v \ is significantly 
lower in C2 compared to Cl [F(l, 8) = 28.501, p < 0.01]. Again, as expected, no significant effects are 
found for K v i for the change in roll or pitch dynamics, or the growth rate. The change in dynamics or the 
growth rate also did not introduce a significant effect in the final visual gain K v 2 . The values for K V 2 are 
significantly lower than K v \ for the conditions with a change in dynamics [F(l, 8) = 69.017, p < 0.001]. In 
these conditions the value for K V 2 more closely resembled the value for K v \ in C2. This was expected for the 
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Figure 11. Means and 95% confidence intervals of the roll pilot model parameters. 
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Figure 12. Means and 95% confidence intervals of the pitch pilot model parameters. 
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conditions with a variation in pitch dynamics, however, not for the conditions with only a variation in roll 
dynamics. This indicates that control behavior in the pitch axis is influenced by the time-varying aircraft 
dynamics in roll. 

Comparing the visual gains between the roll and pitch axes (Fig. 11a with Fig. 12a) it was found that 
the initial visual gain K v i was significantly higher in the pitch axis [F(l, 8) = 7.469, p < 0.05]. The final 
visual gain K v 2 is not significantly different in the roll or pitch axis [F(l, 8) = 0.514, p > 0.05]. 

From Fig. lib it can be observed that in the roll axis the visual lead Tn is significantly higher in C2 
compared to Cl [F(l, 8) = 48.471, p < 0.001]. Based on the characteristics of the aircraft dynamics this was 
an expected result and also supports findings in previous studies. The initial visual lead Tj 1 was not affected 
by a change in dynamics or the growth rate, again an expected result as this parameter describes control 
behavior for the stable dynamics H c 1 in all conditions, except C2. The change in dynamics introduced a 
significant effect in the final lead time constant Tj 2 [F( 2, 16) = 38.755, p < 0.001]. For the conditions with 
a change in roll dynamics, T 12 is higher and resembles T) 1 in condition C2. A comparison of Tn with T 12 
showed a significantly higher value is being found for Tj 2 [F(l, 8) = 345.814, p < 0.001]. 

The effects for the visual lead time constant in the pitch axis were similar to the effects in the roll axis 
(Fig. 12b). In the pitch axis, the lead time constant Tn is also significantly higher for condition C2 compared 
to Cl [F(l, 8) = 190.593, p < 0.001]. The change in dynamics in roll or pitch, or the growth rate did not 
introduce a significant effect in Tn. The change in dynamics did introduce a significant effect in the final 
lead time constant T 12 [F(2, 16) = 50.122, p < 0.001]. The lead time constant for the conditions with a 
change in pitch dynamics approaches the value of the lead time constant found for C2. Comparing Tn with 
Tn in the pitch axis reveals that T) 2 is significantly higher [F(l, 8) = 302.983, p < 0.001]. In the conditions 
with only a change in roll dynamics (C3 and C6), the lead time constant in the pitch axis is significantly 
higher, indicating a cross coupling effect. 

When the lead time constant is compared between the roll and pitch axes (Fig. lib with Fig. 12b), it is 
found that the initial lead time constant Tn is not significantly different between the two axes [F(l, 8) = 
0.884, p > 0.05]. The final lead time constant T) 2 is significantly higher in the pitch axis compared to the 
roll axis [F(l, 8) = 8.720, p < 0.05], 

Figs. 11c and 12c depict the time of maximum growth for the roll and pitch axis, respectively. No 
significant effects were found for the change in roll or pitch dynamics. In both axes the growth rate introduced 
a significant effect, [F(l, 8) = 16.573, p < 0.01] and [F( 1, 8) = 15.529, p < 0.01]. For the higher growth 
rate, representing the step-like change, the time of maximum growth is smaller, indicating that the change 
in pilot control behavior occurs sooner after the change in aircraft dynamics at t = 50 s. No significant 
differences in M were found between the roll and pitch axes [F(l, 8) = 0.099, p > 0.05]. 

The growth rate G of the pilot visual gain and lead time constant is presented in Fig. lid for the roll axis 
and Fig. 12d for the pitch axis. The value of the growth rate had an upper bound of 100 in the estimation 
procedure. In both the roll and pitch axes, the growth rate of the aircraft dynamics did not affect the growth 
rate of the change in pilot dynamics, [F(l, 8) = 1.078, p > 0.05] and [F(l, 8) = 1.344, p > 0.05]. The change 
in dynamics did not introduce a significant effect in the roll axis [F(2, 16) = 3.063, p > 0.05], however in 
the pitch axis the growth rate is significantly higher in the conditions with only a change in roll dynamics 
[F(2, 16) = 60.643, p < 0.001]. Comparing the growth rates between both axes (Fig. lid with Fig. 12d), the 
growth rate in the pitch axis is significantly lower [F(l, 8) = 36.409, p < 0.001]. Note that the error bars 
for the growth rate are relatively large, indicating large variations between subjects. This is most likely an 
effect from the lower accuracy of this parameter in the estimation procedure, as found in Sec. III.B. 

The pilot visual time delay t v is depicted in Figs, lie and 12e for the roll and pitch axes, respectively. 
In the roll axis, the visual time delay is not significantly different between the two conditions with constant 
dynamics (Cl and C2) [F( 1, 8) = 0.259, p > 0.05]. This supports the assumption that this parameter does 
not change over time when the aircraft dynamics change as described in Sec. II. In the pitch axis, the visual 
time delay is significantly different between conditions Cl and C2 [F(l, 8) = 18.950, p < 0.01]. However, 
note that the difference is only 20 ms. In the roll axis, the visual time delay is significantly affected by a 
change in aircraft pitch dynamics [F(l, 8) = 23.923, p < 0.01]. When a change in pitch dynamics occurs, 
the time delay in roll is significantly lower. The change in roll dynamics or the growth rate did not change 
the time delay significantly, [F(l, 8) = 2.685, p > 0.05] and [F(l, 8) = 0.062, p > 0.05]. Comparing the time 
delay between the roll and pitch axes (Fig. lie with Fig. 12e), the time delay in the pitch axis is significantly 
higher [F(l, 8) = 19.753, p < 0.01], 

Figs. Ilf and 12f depict the neuromuscular damping for the roll and pitch axis, respectively. In both 
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axis, the damping is not different between conditions Cl and C2, supporting the assumption to define this 
parameter as a constant parameter, [-F(l, 8) = 1.851, p > 0.05] and [F(l, 8) = 0.229, p > 0.05]. In the roll 
axis, no significant effects were found for the change in roll or pitch dynamics, or the growth rate ([F(l, 8) = 
0.135, p > 0.05], [F(l, 8) = 4.945, p > 0.05] and [F(l, 8) = 0.757, p > 0.05]). The change in dynamics and 
the growth rate also did not introduce a significant effect in the pitch axis ([F(l, 8) = 0.563, p > 0.05], [F(l, 
8) = 1.039, p > 0.05] and [F(l, 8) = 1.060, p > 0.05]). Comparing the neuromuscular damping between the 
two axes, the damping was found to be significantly lower in the pitch axis, [F(l, 8) = 16.467, p < 0.01]. 

The neuromuscular frequency is given in Figs, llg and 12g for the roll and pitch axes, respectively. 
There was no significant difference between conditions Cl and C2, supporting the assumption of a constant 
parameter, [F(l, 8) = 0.706, p > 0.05] and [F(l, 8) = 0.976, p > 0.05]. In the roll axis, the neuromuscular 
frequency was significantly higher when there was a change in roll dynamics present [F(l, 8) = 9.277, p < 
0.05] and significantly lower when there was a change in pitch dynamics present [F(l, 8) = 55.881, p < 
0.001]. The growth rate did not introduce a significant effect, [F(l, 8) = 1.313, p > 0.05]. In the pitch axis, 
no significant effects were found for the change in roll or pitch dynamics, or the growth rate ([F(l, 8) = 
0.053, p > 0.05], [F(l, 8) = 0.226, p > 0.05] and [F(l, 8) = 0.641, p > 0.05]). Comparing the neuromuscular 
frequency between the two axes, the frequency was found to be significantly higher in the pitch axis, [F(l, 
8) = 44.461, p < 0.001]. 

An example of the time-varying pilot visual gain and visual lead time constant as defined by the param- 
eters estimated using MLE is depicted in Fig. 13 for both the roll and pitch axes of condition C8 for one 
pilot. The significant trends identified in Figs. 11 and 12 are clearly visible in the parameter functions. 

(a) Visual gain (b) Visual lead time constant 




Figure 13. Time- varying pilot model parameters (1 pilot, condition C8). 


D. Pilot- Vehicle Open Loop 

The crossover frequency and phase margin of the pilot-vehicle open loop provide insight into the performance 
and stability margins of the pilot-vehicle system. The crossover frequency co c , is the frequency where the 
magnitude of the open- loop response has a magnitude of 1. The corresponding phase margin ip is the phase 
difference with -180 deg at the crossover frequency. Fig. 14 provides an example of the open-loop frequency 
response for the roll axis of condition C8 for one pilot. A clear increase in open-loop gain and decrease 
in effective time-delay r e can be seen after the change in aircraft dynamics. The result is an increase in 
crossover frequency and a decrease in phase margin, as depicted in Fig. 15. This figure also provides the 
time-varying crossover frequency and phase margin for the pitch axis. 

VI. Discussion 

A parameter estimation technique based on MLE was used to identify a pilot model with time-dependent 
parameter functions characterizing time-varying pilot control behavior. The technique was used to estimate 
the pilot model parameters using experimental data from a multi-axis disturbance-rejection task. 

Several assumptions were made with respect to the time-varying pilot model. The first was the type 
of parameter function to describe the time-varying parameters, which was assumed to be equivalent to the 
parameter function used to vary the aircraft dynamics over time. Given the aircraft dynamics used in the 
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Figure 14. Time- varying open-loop frequency response for the roll axis (1 pilot, condition C8). 
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Figure 15. Crossover frequency and phase margin as a function of time (1 pilot, condition C8). 


experiment, this was a straightforward choice. However, if time-varying pilot control behavior needs to be 
identified when little is known about the characteristics of the time- varying controlled dynamics, the selection 
of the most appropriate parameter function can be more difficult. It requires incorporating the available 
knowledge of the time- varying aspect of the controlled dynamics and careful inspection of the measured data. 

The second assumption is the introduction of an equivalent time of maximum growth and growth rate 
for both parameter functions in the pilot model. This means that the time-varying parameters of the model 
change at the same time and with the same rate with only their initial and final values differing. This reduces 
the number of parameters to be estimated, resulting in more accurate parameter estimates. However, this 
assumption may not hold for all applications and must be considered in conjunction with the selection of 
the parameter function type. 

The final assumption is to only let the pilot equalization vary over time, in this case the visual gain and 
visual lead time constant. The remaining parameters - the visual time delay, the neuromuscular damping, and 
the neuromuscular frequency - were assumed to be constant over time. This assumption was also necessary to 
reduce the number of parameters to be estimated in the MLE procedure, increasing the accuracy. Previous 
research has shown that pilots mainly adapt their equalization dynamics when controlling the different 
controlled dynamics used in the current study. 15 The assumption was validated by comparing constant pilot 
control behavior between two reference conditions with constant aircraft dynamics, Cl and C2. The aircraft 
dynamics in these conditions are the same as the initial and final aircraft dynamics in conditions C3-C8. 
Comparing the estimated pilot model parameters between Cl and C2 shows that only the visual gain and 
lead time constant were significantly different between the two conditions. 

To validate the accuracy of the pilot model, the VAF was calculated for the roll and pitch control signal. 
The VAF was generally around 70 % for the roll axis and around 82 % for the pitch axis. These values are 
slightly lower compared to values found in single-axis experiments. This was an expected result, since in 
multi-axis control tasks the pilot’s attention needs to be divided between the two axes, resulting in higher 
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levels of pilot remnant. The higher VAF for the pitch axis indicates that pilot control behavior was more 
linear in the pitch axis (that is, remnant levels were lower) compared to the roll axis. This could be due to 
the fact that a pitch disturbance is more visible than a roll disturbance due to the nature of the primary 
flight display used in the experiment (a roll disturbance is depicted by an angle deviation, while a pitch 
disturbance is depicted by a vertical displacement), or simply by the fact that pilots paid more attention to 
pitch compared to roll. 

Furthermore, as hypothesized in Sec. IV. B, the pilot model VAF for the unstable double-integrator-like 
dynamics H c 2 was higher compared to the pilot model VAF for the single-integrator-like dynamics H c 1 . This 
was also found in previous experiments. 17 Finally, for the conditions with time- varying aircraft dynamics 
(C3-C8), a constant pilot model with 5 parameters H p (s) was identified in addition to the time- varying pilot 
model with 9 parameters H p (s,t), to verify if the added time-varying parameter functions provide a higher 
overall model accuracy. This was indeed the case, for all conditions with time-varying aircraft dynamics, the 
time-varying pilot model had a significant higher VAF compared to the constant model. The difference in 
VAF was greater in the pitch axis (approximately 7 %) compared to the roll axis (approximately 3 %). 

The accuracy of the MLE procedure and the selected pilot model structure was also analyzed in a Monte 
Carlo analysis before the start of the experiment. This analysis showed that the final values of the parameter 
functions and the growth rate have the highest bias and variance, and that bias and variance increase when 
pilot remnant intensity increases. The lower accuracy of the growth factor can be explained by the fact that 
this value is only affected by data from a small time window. In addition, the growth rate does not have an 
upper limit; that is, for values of approximately 10 to infinity, the parameter function is a step function. For 
this reason, the growth rate had an upper bound of 100 in the estimation procedure. 

The independent variables - the changing aircraft dynamics in the roll and pitch axes, and the growth 
rate of the dynamics - introduced significant adaptations in pilot control behavior. Mainly the parameters 
defining the time-dependent sigmoid functions were significantly affected. The main effect was in the axis of 
the changing aircraft dynamics. However, significant effects were also found in the axis other than the axis 
of the changing aircraft dynamics. This indicates some cross coupling in pilot control behavior between the 
two axes, and might warrant the addition of a cross coupling term in a pilot model predicting multi-axis 
control behavior, even when the controlled dynamics in the two axis are not coupled. 

When comparing the results across conditions between the roll and pitch axes, it was found that perfor- 
mance is higher and control activity is lower in the pitch axis. In addition, the average initial visual gain, 
final visual lead time constant, time delay, and neuromuscular frequency were all significantly higher for 
the pitch axis, while the neuromuscular damping was significantly lower. Note that the aircraft dynamics 
and forcing functions are equivalent in both controlled axes, indicating that these effects were caused by 
differences in display and joystick characteristics between the two axis, or differences in human perception 
and control processes between the two axes. 

Most studies on pilot control behavior use control tasks in a single axis and assume constant pilot control 
behavior. These assumptions reduce the applicability of the results from these studies to real piloting tasks, 
as these are often performed in multiple axes and are not constant over time. In this paper time-varying 
pilot control behavior was accurately identified using data from a multi-axis control task, thus expanding the 
scope of pilot control behavior research. Further research is required to investigate the use of different types 
of parameter functions to allow for more complex parameter variations in time. In addition, time-variations 
in parameters other than the pilot equalization, such as time-varying perceptual time delays, should be 
investigated. 


VII. Conclusions 

A maximum likelihood parameter estimation method was used to identify a pilot model with time- 
dependent sigmoid functions to characterize time-varying human control behavior. An experiment was 
performed with 9 general aviation pilots performing a simultaneous roll and pitch control task with time- 
varying aircraft dynamics. In 8 different conditions the controlled axis containing the time- varying dynamics 
and the growth factor of the change in dynamics was varied. The MLE method used was able to estimate 
most parameters of the time- varying pilot model from the experiment data accurately. However, the growth 
factor of the sigmoid functions had the lowest estimation accuracy, as it is unbounded for higher values due 
to the nature of the sigmoid function. Time- varying pilot control behavior in both the roll and pitch axes was 
significantly affected by the axis containing the time-varying aircraft dynamics and the growth factor. The 
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main effect was found in the axis that contained the time-varying aircraft dynamics. However, significant 
effects were also found in the opposite axis. This indicates that cross coupling exists in pilot control behavior 
between the two axes. Although the same aircraft dynamics and forcing functions are used in both axes of 
control, pilot performance was higher and control activity was lower in pitch compared to roll. This was also 
reflected by significantly different pilot control behavior in both axes. 
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